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Abstract. We propose a computational method to simulate anomalous self-diffusion 
in a simple liquid. The method is based on a molecular dynamics simulation on 
which we impose the following two conditions: firstly, the inter-particle interaction is 
described by a soft-core potential and secondly, the system is forced out of equilibrium. 
The latter can be achieved by subjecting the system to changes in the length scale 
at intermittent times. In many respects, our simulation system bears resemblance 
to slowly driven sandpile models displaying self-organised criticality. We find non- 
Gaussian single time step displacement distributions during the out-of-equilibrium time 
periods of the simulation. 
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1. Introduction 

Ever since Einstein made geometric Brownian motion famous [Tj, a simple liquid has 
been the prime example of a system in which normal diffusion occurs. Gaussian 
statistics, which characterizes normal diffusion, is extensively utilized in various fields 
of research. Anomalous, non-Gaussian diffusion is characterized by a probability 
distribution function with power-law tails P(z > z ) ~ \z\~ a . The limiting distribution 
of the sum of independent variables with this property is a Levy distribution [21 |3l H] . 
It is well known that anomalous diffusion, in the general sense, can occur in both 
equilibrium and non-equilibrium contexts [SJ E] ■ Amongst the fields in which examples 
of non- Gaussian dynamics are found are economics [3 El [9], biology [TUJ [TTJ and solid 
state physics [P2] . 

The above observations have motivated extensive efforts to investigate the 
underlying mechanisms of non-Gaussian dynamics. There are numerous methods 
to achieve anomalous diffusion either in lower-dimensional systems [13J, in media 
P~H [15], in networks [16] or in turbulence [T7]. To our knowledge, however, in 
literature there is no mention of a simple computational method to achieve conditions 
of anomalous self-diffusion in a 3D liquid. Here, we wish to propose such a method. 
We present a robust technique to alter the dynamics of a molecular dynamics (MD) 
simulation of a simple liquid in such a manner that anomalous self-diffusion emerges. 
These anomalous properties emerge in the average one- dimensional single time step 
displacement distributions P(Ax). In particular, we propose to combine a MD 
simulation technique with a soft-core inter-particle potential and out-of-equilibrium 
conditions. 

MD is a very powerful simulation technique that is based on (numerically) solving 
the equations of motion for many interacting particles at carefully chosen regular time 
intervals. As a consequence, a MD simulation allows one to investigate many aspects of 
liquids and other systems such as self-diffusion, phase diagrams, absorption of particles 
and viscosity [T51 [T51 I2"0] |2"T] . It is well known that in classical MD simulations with 
a Lennard- Jones (LJ) inter-particle potential, the mean-square displacement (MSD) 
(Ar 2 (t)) has a linear time dependence and that P(Ax) is a Gaussian distribution. 

2. Formalism 

When attempting to create conditions under which anomalous self-diffusion emerges, the 
hard core of the LJ potential, Ulj(t) = e ((ro/r) 12 — 2 (r /r) 6 ), poses a real challenge. 
Indeed, in a liquid with conditions of anomalous self-diffusion, P(Ax) obtains heavy 
tails and the particles can occasionally traverse a relatively large distance during a 
single time step. In the MD simulations, those particles will momentarily travel with 
an exceptionally large speed. Due to the finite time resolution of the simulation, they 
can penetrate the hard core of the LJ potential and attain a velocity that is much larger 
than the average velocities in the simulation. This, in turn, dramatically increases the 
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probability of particles to enter the hard core of other particles and ignites a chain 
reaction that makes the simulation to go out of control. In order to remedy this, we 
have renormalized the short-range part of the Lennard- Jones potential and introduced 
the following soft-core potential [22] : 



U SC {r) = — -r~, - U A exp 

1 + exp A (r — R R ) 



25\ 



(1) 



The parameters A, R R , Ra, U a, 5a were optimized so as to match the medium- and long- 
range part of ULj{r). The parameter Ra determines the length scale in the system and 
we take Ra = r = 3.405 x 10 _10 m, the diameter of Argon in a LJ potential. The 
adopted units for energy and mass are also determined for a simulation that involves 
Argon: m = 6.63 x 10~ 26 kg and E = 1.65 x 10~ 21 J. The parameter H determines 
the hardness of the core of the potential. The softer the core, the more penetrable the 
particles become. 

First, we wish to discriminate between the dynamical properties of the mono-atomic 
liquid for Uij(r) and Usc( r )- To this end, the potential Usc( r ) of Q is used in a classical 
MD simulation program. We investigate the self-diffusion properties, the velocity auto- 
correlation function (VACF) and the radial distribution function (RDF). The VACF 
tracks the persistence of the motion over time and the RDF provides information about 
the relative particle positions, providing a means of discriminating between the liquid, 
gaseous or solid phase. The thermodynamic properties like temperature, potential and 
kinetic energy, are also monitored. In order to solve the equations of motion, we adopt 
the Verlet algorithm, which is a symplectic integrator. 




Figure 1: The MSD as a function of simulation time for different soft-core parameters 
H and for a LJ potential. The simulation has an initial density of 0.5 and an initial 
temperature of 0.7. The density of the system is given in particles per unit of volume 
and the temperature is given in system units. All distances are expressed in system 
units, which are determined by Ra = r Q . The simulations are performed with 8788 
particles, and the time step is determined by aoo !^= 1 . 
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Figure [T] shows the time dependence of the MSD for various values of H of the 
soft-core potential and for a LJ potential. For large values of H, the soft-core potential 
of ([TJ is qualitatively similar to Ulj(t). For all values of if, the dynamics of the system 
lead to normal diffusion ((Ar 2 ) ~ t). The simulation results indicate that the diffusion 
coefficient D ((Ar 2 ) = 6Dt) depends on the softness of the potential. The harder the 
short-range part the lower D will be, as the collisions approach a hard-sphere interaction. 
The qualitative features of the VACF, RDF, energy and temperature obtained with 
Usc{ r ) resemble those of a simulation with a LJ potential. There are some differences 
in the computed observables and the RDF, for example, reflects that the soft core allows 
penetration. The positional structure, typical for a liquid, is still visible in the RDF. The 
most important aspect, however, is the ubiquitous presence of the Gaussian statistics 
in the self-diffusion properties. 



3. Out-of-equilibrium simulation 

We now turn to simulations under non-equilibrium conditions. Under equilibrium 
conditions, the energy is conserved and the temperature fluctuates mildly around a 
certain value. We introduce non-equilibrium conditions by driving the system and 
modifying the inter-particle interaction. This can be achieved by rescaling the radial 
distances in the soft-core interaction Uscij") —> Usc(^ r ) with A < 1. This is equivalent 
to an effective increase of the size of the molecules. This allows the dynamics of the 
system to be changed dramatically, because particles that were attracting each other 
end up repelling each other due to the driven change in the inter-particle interaction 
range. As the imposed changes in the inter-particle interactions occur under conditions 
of constant density, the system develops regions of high energy density. Through the 
dynamics of the system, local energy surplus dissipates into sizeable kinetic energy and 
an increase in the temperature of the system is observed. Figure [2] shows the evolution 
of the temperature with time for a simulation that undergoes a rescaling of the soft- 
core interaction of the particles. It is clear that the driven change in the radius of the 
molecules has a large effect on the temperature. After rescaling the radial distances, 
it takes of the order of a few hundred time steps for the temperature to reach a new 
equilibrium value. 

To illustrate the effect of the radial rescalings on the internal dynamics of the 
system, the potential energy fluctuation of the system during one simulation step is 
shown in figure [3j We consider results for Ulj(t) and Usc{ r ) under typical equilibrium 
conditions and a situation just after a radial rescaling. For every particle the difference 
in potential energy with respect to the previous configuration is shown 

AE pot ( ri ) =Y i U (r - (t + At)) - U (r v (t)) . (2) 
The panels of figure [3] represent a projection of a slice (Vz : \Zi\ < 0.5) onto the 



xy-plane. Figure 3a indicates that through a sudden rescaling of the radial distances 



one creates regions in which the total amount of potential energy gain is much larger 
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Figure 2: Temperature as a function of time for a simulation in which a radial rescaling 
with A = 0.75 takes place every 500 time steps. When the temperature reaches 50, the 
radii and the temperature are reset to their original value (T init = 0.7). 




Figure 3: A typical spatial distribution of the potential energy changes AE pot {vi) in one 
time step. We show the projection onto the xy-plane for \zi\ < 0.5 under conditions of 
(a) a soft-core potential just after rescaling with A = 0.7 (neq), (b) a LJ potential just 
after rescaling with A = 0.7 (neq), (c) a soft-core potential during equilibrium (eq), (d) 
a LJ potential during equilibrium (eq). 



than the average value. Figure 3b shows that the hard core of Ulj(t) results in values 
of AEp t as high as 800, whereas this is not seen in figure 3a With energy fluctuations 
of this size, the velocities of the particles attain values that are not compatible with a 
finite time-step. Figure 3c and figure 3d show a similar type of projection for typical 
equilibrium conditions. The scale of AE pot is clearly much smaller than in the non- 



equilibrium situations of figures 3a and 3b 



The regular rescaling of the inter-particle distances drives the system's 
thermodynamic properties such as the temperature and the energy away from 
equilibrium, i.e. mild fluctuations around a constant value. As a result, the simulation 
resembles conditions encountered in systems which display self-organized criticality 
(SOC) [23]. These systems are characterized by a driving and a relaxation mechanism. 
In our studies, the increase of the radius is the driving factor that injects potential 
energy at various positions in the system and the conversion from potential to kinetic 
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energy is the relaxation mechanism. The balance between the injected and dissipated 
energy is linked through the local Newtonian dynamics that conserves total energy 
Another characteristic feature of SOC is that the system must be driven slowly. For 
our purposes, this translates into changing A after sufficiently long time intervals. For 
a LJ-potential, the driven changes in the distance scale of the inter-particle interaction 



amount to very large potential energy changes (Figure 3b ) that are not compatible with 
slowly driving the system. 

As repeatedly increasing the radii is not an attractive option (because of the finite 
size of the simulation system), after some time we reset the system's original temperature 
and particle radius. If the temperature exceeds a certain threshold, the radii are rescaled 
to their original value and the velocities are rescaled so that the starting temperature is 
restored, see figure [2] 

4. Results 

We identify the time intervals with a varying temperature in figure [2] as non-equilibrium 
conditions. Now, we study the self-diffusion properties of the liquid under those non- 
equilibrium conditions. For a single time-step, the mean, standard deviation and 
kurtosis of P(Ax), P(Ay) and P(Az) are calculated. The standard deviation (cr) and 
kurtosis (k) are defined in the standard fashion, 



a(t) 



1 N 



£(A^) - Ax(tW (3) 



with fi4 the fourth moment about the mean. This definition ensures that for a Gaussian 
distribution, the kurtosis vanishes. 

To normalise the width of the different distributions, the displacements are divided 
by the standard deviation of the distribution at some time instance. This normalisation 



enables one to compare the P(Ax/o~) that are obtained at different times. Figure 4a 
compares P(Ax/a) for a typical equilibrium time instance with one obtained at a 
representative non-equilibrium time instance. Under non-equilibrium conditions, the 
tails of P(Ax/o~) are considerably fatter than under equilibrium conditions. Moreover, 
a higher concentration of particles in the centre of the distribution is observed. These 
characteristics are typical for a leptokurtic distribution. The kurtosis of the out-of- 
equilibrium P(Ax/a) are larger than the typical noise level under equilibrium conditions. 

We have found that after effectively enlarging the particles, anomalous 
characteristics can be found in the single time step displacement distributions P(Ax/o~). 
After rescaling the interaction four times with A = 0.75 and 500 time steps between the 
different rescalings (sufficient to reach equilibrium again), we have obtained anomalous 
distributions in P{Ax/a). To establish the non-Gaussian shape of the tail of P(Ax/a) 
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Figure 4: P(Ax/o~) for two typical equilibrium (Gaussian) and non-equilibrium 
(anomalous) situations with initial density p = 0.5, temperature T = 0.7. The upper 
panels are for A = 0.7 after four rescalings. The lower panels are for A = 0.75 after 
four rescalings. The right panels are a fit of P(Ax/a > 1.5) with Q. The best fit 
parameters are a = 792 and b = 0.14 for the upper right panel and a = 628 and b = 0.13 
for the lower right panel. 



under non-equilibrium conditions, we have fitted it with an exponential: 



,Ax. 



a 



aexp 



Ax 



a 



Ax 



> 1.5 



a 



(5) 



As can be seen in figure 4b and 4d, the distributions are well fitted by ([5]). This result 
confirms that the P(Ax/a) have fatter tails than a Gaussian distribution. 

The distributions of figures 4a and 4c are obtained by taking data during one time 
step. Summing these distributions for different time steps results in better statistics. 
Figure 5a shows the result for P(Ax/a) after summation of 50 distributions during 
anomalous and normal (Gaussian) simulation conditions. Remark that > 4 events 
are one order of magnitude more likely under anomalous (non-equilibrium) conditions 
than under Gaussian (equilibrium) conditions. We find a fair amount of 5a events 
and some rare 8a events. The normalized cumulative distribution of figure [5] clearly 
illustrates that the self-diffusion properties of the liquid are distinctive during the non- 
equilibrium and equilibrium periods of the simulation. In non-equilibrium conditions, 
small Ax I a are more likely, medium Ax /a less likely and large Ax /a far more likely. 
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Figure 5: Left: Summation of 50 'anomalous' and 'Gaussian' distributions of P(Ax/a). 
Right: Normalized cumulative distribution of the sum of 50 'anomalous' and 'Gaussian' 
distributions P(Ax/a). 



We now wish to study the trajectories of the individual particles during the complete 
simulation that alternates equilibrium with non-equilibrium conditions. To this end, we 
selected two particles: particle #1524 for which \Ax/a\ does not exceed 5 during the 
simulation (since this limit isn't attainable in a Gaussian simulation regime) and particle 
#3329 for which this is not the case. Figure [6] illustrates that the random character of 
the trajectories for both particles is maintained over the simulation time. 



2000 




-1500 -1000 -500 500 1000 

I Ax/a 

Figure 6: Projection on the xy-p\a.iae of J2t Ar(t)/er(t) °f particles #1524 and #3329 
during a simulation of 100000 time-steps. 

Figure [7] shows Ax /a as a function of time for these particles. The "anomalous" 
behaviour of particle #3329 is confined to the time period 32000 - 33000. By contrast, 
no anomalous behaviour is discernible in the behaviour of particle #1524. 
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Figure 7: Single time step displacements (Ax(t)/a(t)) of particles #1524 and #3329 
during a simulation of 100000 time steps. 

To establish the robustness of our technique to generate conditions of anomalous 
self-diffusion in a mono-atomic liquid, we have investigated its dependence on the 
parameters of the simulations. The non-equilibrium conditions are determined by the 
size of the rescaling parameter A and the time intervals between two subsequent radial 
rescalings. During a simulation of 100000 time-steps, the system goes through different 
periods of non-equilibrium conditions. The kurtosis of P(Ax/cr) is computed for every 
time-step and the maximal kurtosis is saved for every set of parameters. In figure [8] the 
maximum kurtosis of the simulation is plotted as a function of \ and the time interval 
between two rescalings. The time intervals are chosen between 100 and 1500 time-steps, 
since equilibrium is always reached after 1500 time-steps. A is confined to 1 < 1/A < 1.4, 
which means that the density change in one step is limited to 1 < Ap < 2.7. Figure [8] 
clearly shows that the anomalous character of P(Ax/o~) remains present independent 
of the rescaling parameters. The time interval between two subsequent rescalings (r) 
has only a minor influence on the maximal kurtosis (k max ). The adopted value of A, on 
the other hand, has a larger influence on k max , but the anomalous characteristics are 
present in the entire A interval. 

The same robustness applies to variations in the initial density and temperature of 
the system, provided that they generate a system in the liquid phase. The short range 
order and particle mobility typical for a liquid are necessary conditions. They generate 
the required balance between mobility and interaction, allowing the system to dissipate 
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Figure 8: Maximal kurtosis (k max ) of a simulation as a function of the rescaling 
parameter \ and the time interval r between two subsequent rescalings. 

The robustness of the anomalous character of the self-diffusion properties under 
non-equilibrium conditions is a very useful result. It indicates that the qualitative 
features of the self-diffusion properties of the system are rather insensitive to the 
two parameters (A and r) that characterise the non-equilibrium behaviour. As a 
consequence, we can ascertain that it is the internal dynamics of the system that causes 
the non-Gaussian properties of P(Ax/a). 

We have also tested the robustness of our results to changes in the potential. To this 
end, we have performed simulations with the following potential: a LJ for ry > 0.9 and 
a polynomial ar 6 + b for < 0.9 [24]. This potential has the long-range LJ properties 
and a soft core. We have obtained non-Gaussian distributions almost identical to those 
obtained with Usc( r )- 

A further test of the robustness of our technique can be done in dimensions other 
than three. In two dimensions we obtain inherent anomalous self-diffusion of the system, 
as expected [3 [13]. Figure [9] shows the result of a simulation in four dimensions. We have 
used the same technique as the one adopted for the 3D results of Figure [5| From figure [9] 
it is clear that during the non-equilibrium time periods the self- diffusive properties 
are non-Gaussian. This provides further evidence for the robustness of the proposed 
technique. 

5. Conclusion 

In summary, we have presented a computational method, based on out-of-equilibrium 
MD simulations with a soft-core potential, that generates dynamical conditions of 
anomalous diffusion for the distribution of the one-dimensional displacements of the 
particles. This behaviour arises because the driving mechanism, i.e. the potential 
energy artificially injected into the system by increasing the radius of the particles, 
generates regions of increased potential energy. The dissipation of this potential energy 
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Figure 9: P(Ax/a) for a four-dimensional simulation with 4375 particles, summed 
over 350 time instances, for two typical equilibrium (Gaussian) and non-equilibrium 
(anomalous) situations, with A = 0.8 and r = 1000. 



under conditions of constant energy results in a small amount of very fast particles. For 
a broad range of the parameters involved, our technique generates anomalous diffusion 
in the simulation. In this way we have created a simple and elegant method to achieve 
anomalous self-diffusion in an interacting system. This emergence of global behaviour 
that cannot be determined from local properties is also a property of self-organized 
criticality, to which our model bears similarities. 
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